#***************************************************************#
#                                                               #
#                                                               #
#                 REPLICATION FILES FOR:                        #
#                                                               #
#                  Endogenous Democracy:                        #
#  Causal Evidence from the Adoption of Potato in the Old World #                           
#                                                               #
#                                                               #          
#                                                               #
#                                                               #
#***************************************************************#


#***************#
# Load Packages #
#***************#

library (systemfit)
library (AER)
library (foreign)
library (plm)
library (survival)
library (apsrtable)
library (DataCombine)
library(ivpack)
library(utils)
library(ggplot2)

#***********************#
# Set Working Directory #
#***********************#

Path <- ""

#setwd(Path)

#***********************#
#       Load data       #
#***********************#

Potato <- read.table("potato_data.txt", sep="\t", header=T)


#***********************#
#        Table 1        #
#***********************#

#**********#
# Column 1 #
#**********#

summary(t1_column1 <- ivreg(lead.polity ~ city_pop_share + as.factor(country) + as.factor(year) - 1 | . - city_pop_share + ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(t1_column1, t1_column1$model$`as.factor(country)`)

summary(first1_column1 <- lm(city_pop_share ~ ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900 + as.factor(country) + as.factor(year), data=Potato), diagnostics=TRUE)

cluster.robust.se(first1_column1, first1_column1$model$`as.factor(country)`)


#**********#
# Column 2 #
#**********#

summary(t1_column2 <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1100 + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + ln_oworld_1850 + ln_oworld_1900 + as.factor(country) + as.factor(year) - 1 | . - city_pop_share + ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(t1_column2, t1_column2$model$`as.factor(country)`)

summary(first1_column2 <- lm(city_pop_share ~ ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900 + + ln_oworld_1100 + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + ln_oworld_1850 + ln_oworld_1900 + as.factor(country) + as.factor(year), data=Potato), diagnostics=TRUE)

cluster.robust.se(first1_column2, first1_column2$model$`as.factor(country)`)


#**********#
# Column 3 #
#**********#

summary(t1_column3 <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1100 + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + ln_oworld_1850 + ln_oworld_1900 + ln_elevation_1100 + ln_elevation_1200 + ln_elevation_1300 + ln_elevation_1400 + ln_elevation_1500 +  ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + ln_tropical_1100 + ln_tropical_1200 + ln_tropical_1300 + ln_tropical_1400 + ln_tropical_1500 + ln_tropical_1600 + ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + ln_tropical_1850 + ln_tropical_1900 + ln_rugged_1100 + ln_rugged_1200 + ln_rugged_1100 + ln_rugged_1200 + ln_rugged_1300 + ln_rugged_1400 + ln_rugged_1500 + ln_rugged_1600 + ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + ln_rugged_1900 + as.factor(country) + as.factor(year) - 1 | . - city_pop_share + ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(t1_column3, t1_column3$model$`as.factor(country)`)

summary(first1_column3 <- lm(city_pop_share ~  + ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900 + ln_oworld_1100 + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + ln_oworld_1850 + ln_oworld_1900 + ln_elevation_1100 + ln_elevation_1200 + ln_elevation_1300 + ln_elevation_1400 + ln_elevation_1500 + ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + ln_tropical_1100 + ln_tropical_1200 + ln_tropical_1300 + ln_tropical_1400 + ln_tropical_1500 + ln_tropical_1600 + ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + ln_tropical_1850 + ln_tropical_1900 + ln_rugged_1100 + ln_rugged_1200 + ln_rugged_1100 + ln_rugged_1200 + ln_rugged_1300 + ln_rugged_1400 + ln_rugged_1500 + ln_rugged_1600 + ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + ln_rugged_1900 + as.factor(country) + as.factor(year), data=Potato), diagnostics=TRUE)

cluster.robust.se(first1_column3, first1_column3$model$`as.factor(country)`)


#*********************************#
#                                 #
#        Online Appendices        #
#                                 #
#                                 #
#*********************************#

#***********************#
#       Appendix D:     #
#     Robustness Checks #
#******************************#
# Restricting pre-1800 periods #
#******************************#

#****Models from Column 1******#

#**Including all observations**#

summary(c1_noexcl <- ivreg(lead.polity ~ city_pop_share
                           + as.factor(country) + as.factor(year) - 1
                           | . - city_pop_share + ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300
                           + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                             ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_noexcl, c1_noexcl$model$`as.factor(country)`)


#**Including observations posterior to 1100**#

summary(c1_excl1100 <- ivreg(lead.polity ~ city_pop_share + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1200 + ln_wpot_1300
                             + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_excl1100, c1_excl1100$model$`as.factor(country)`)

#**Including observations posterior to 1200**#

summary(c1_excl1200 <- ivreg(lead.polity ~ city_pop_share
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1300
                             + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_excl1200, c1_excl1200$model$`as.factor(country)`)

#**Including observations posterior to 1300**#

summary(c1_excl1300 <- ivreg(lead.polity ~ city_pop_share
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_excl1300, c1_excl1300$model$`as.factor(country)`)

#**Including observations posterior to 1400**#

summary(c1_excl1400 <- ivreg(lead.polity ~ city_pop_share + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_excl1400, c1_excl1400$model$`as.factor(country)`)

#**Including observations posterior to 1500**#

summary(c1_excl1500 <- ivreg(lead.polity ~ city_pop_share 
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_excl1500, c1_excl1500$model$`as.factor(country)`)

#**Including observations posterior to 1600**#

summary(c1_excl1600 <- ivreg(lead.polity ~ city_pop_share 
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_excl1600, c1_excl1600$model$`as.factor(country)`)

#**Including observations posterior to 1700**#

summary(c1_excl1700 <- ivreg(lead.polity ~ city_pop_share 
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

cluster.robust.se(c1_excl1700, c1_excl1700$model$`as.factor(country)`)

#**Including observations posterior to 1800**#

summary(c1_excl1750 <- ivreg(lead.polity ~ city_pop_share +
                               as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1750 + ln_wpot_1800
                             + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Final results with corrected SE**#

c1_noexcl <- cluster.robust.se(c1_noexcl, c1_noexcl$model$`as.factor(country)`)
c1_excl1100 <- cluster.robust.se(c1_excl1100, c1_excl1100$model$`as.factor(country)`)
c1_excl1200 <- cluster.robust.se(c1_excl1200, c1_excl1200$model$`as.factor(country)`)
c1_excl1300 <- cluster.robust.se(c1_excl1300, c1_excl1300$model$`as.factor(country)`)
c1_excl1400 <- cluster.robust.se(c1_excl1400, c1_excl1400$model$`as.factor(country)`)
c1_excl1500 <- cluster.robust.se(c1_excl1500, c1_excl1500$model$`as.factor(country)`)
c1_excl1600 <- cluster.robust.se(c1_excl1600, c1_excl1600$model$`as.factor(country)`)
c1_excl1700 <- cluster.robust.se(c1_excl1700, c1_excl1700$model$`as.factor(country)`)
c1_excl1750 <- cluster.robust.se(c1_excl1750, c1_excl1750$model$`as.factor(country)`)

#****Models from Column 2******#

#**Including all observations**#

summary(c2_noexcl <- ivreg(lead.polity ~ city_pop_share  + ln_oworld_1100 + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + 
                             ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                             ln_oworld_1850 + ln_oworld_1900
                           + as.factor(country) + as.factor(year) - 1
                           | . - city_pop_share + ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1100**#

summary(c2_excl1100 <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1200 + ln_wpot_1300
                             + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1200**#

summary(c2_excl1200 <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1300
                             + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1300**#

summary(c2_excl1300 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1400 + ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1400**#

summary(c2_excl1400 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1500**#

summary(c2_excl1500 <- ivreg(lead.polity ~ city_pop_share + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share +
                               ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1600**#

summary(c2_excl1600 <- ivreg(lead.polity ~ city_pop_share  + 
                               ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1700**#

summary(c2_excl1700 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1750**#

summary(c2_excl1750 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1800 + 
                               ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Final results with corrected SE**#

c2_noexcl <- cluster.robust.se(c2_noexcl, c2_noexcl$model$`as.factor(country)`)
c2_excl1100 <- cluster.robust.se(c2_excl1100, c2_excl1100$model$`as.factor(country)`)
c2_excl1200 <- cluster.robust.se(c2_excl1200, c2_excl1200$model$`as.factor(country)`)
c2_excl1300 <- cluster.robust.se(c2_excl1300, c2_excl1300$model$`as.factor(country)`)
c2_excl1400 <- cluster.robust.se(c2_excl1400, c2_excl1400$model$`as.factor(country)`)
c2_excl1500 <- cluster.robust.se(c2_excl1500, c2_excl1500$model$`as.factor(country)`)
c2_excl1600 <- cluster.robust.se(c2_excl1600, c2_excl1600$model$`as.factor(country)`)
c2_excl1700 <- cluster.robust.se(c2_excl1700, c2_excl1700$model$`as.factor(country)`)
c2_excl1750 <- cluster.robust.se(c2_excl1750, c2_excl1750$model$`as.factor(country)`)


#****Models from Column 3******#

#**Including all observations**#

summary(c3_noexcl <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1100 + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + 
                             ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                             ln_oworld_1850 + ln_oworld_1900 + ln_elevation_1100 + ln_elevation_1200 + ln_elevation_1300 + ln_elevation_1400 + ln_elevation_1500 + 
                             ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + 
                             ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + ln_tropical_1100 + ln_tropical_1200 + ln_tropical_1300 + 
                             ln_tropical_1400 + ln_tropical_1500 + ln_tropical_1600 + 
                             ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + 
                             ln_tropical_1850 + ln_tropical_1900 + 
                             ln_rugged_1100 + ln_rugged_1200 + ln_rugged_1300 + ln_rugged_1400 + ln_rugged_1500 + ln_rugged_1600 + 
                             ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                             ln_rugged_1900
                           + as.factor(country) + as.factor(year) - 1
                           | . - city_pop_share + ln_wpot_1100 + ln_wpot_1200 + ln_wpot_1300
                           + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                             ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1100**#

summary(c3_excl1100 <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1200 + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + ln_elevation_1200 + ln_elevation_1300 + ln_elevation_1400 + ln_elevation_1500 + 
                               ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + ln_tropical_1200 + ln_tropical_1300 + 
                               ln_tropical_1400 + ln_tropical_1500 + ln_tropical_1600 + 
                               ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + 
                               ln_tropical_1850 + ln_tropical_1900 + 
                               ln_rugged_1200 + ln_rugged_1300 + ln_rugged_1400 + ln_rugged_1500 + ln_rugged_1600 + 
                               ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1200 + ln_wpot_1300
                             + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1200**#

summary(c3_excl1200 <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1300 + ln_oworld_1400 + ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + ln_elevation_1300 + ln_elevation_1400 + ln_elevation_1500 + 
                               ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + ln_tropical_1300 + 
                               ln_tropical_1400 + ln_tropical_1500 + ln_tropical_1600 + 
                               ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + 
                               ln_tropical_1850 + ln_tropical_1900 + 
                               + ln_rugged_1300 + ln_rugged_1400 + ln_rugged_1500 + ln_rugged_1600 + 
                               ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1300
                             + ln_wpot_1300 + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1300**#

summary(c3_excl1300 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1400 + ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + ln_elevation_1400 + ln_elevation_1500 + 
                               ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + 
                               ln_tropical_1400 + ln_tropical_1500 + ln_tropical_1600 + 
                               ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + 
                               ln_tropical_1850 + ln_tropical_1900 + 
                               ln_rugged_1400 + ln_rugged_1500 + ln_rugged_1600 + 
                               ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1400 + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1400**#

summary(c3_excl1400 <- ivreg(lead.polity ~ city_pop_share + ln_oworld_1500 + 
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + ln_elevation_1500 + 
                               ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + 
                               ln_tropical_1500 + ln_tropical_1600 + ln_tropical_1700 + ln_tropical_1750 + 
                               ln_tropical_1800 + ln_tropical_1850 + ln_tropical_1900 + 
                               ln_rugged_1500 + ln_rugged_1600 + 
                               ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1500 + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1500**#

summary(c3_excl1500 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1600 + ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + 
                               ln_elevation_1600 + ln_elevation_1700 + ln_elevation_1750 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + ln_tropical_1600 + 
                               ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + 
                               ln_tropical_1850 + ln_tropical_1900 + 
                               ln_rugged_1600 + 
                               ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1600 + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1600**#

summary(c3_excl1600 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1700 + ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + 
                               ln_elevation_1700 + ln_elevation_1750 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + 
                               ln_tropical_1700 + ln_tropical_1750 + ln_tropical_1800 + 
                               ln_tropical_1850 + ln_tropical_1900 + 
                               ln_rugged_1700 + ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share + ln_wpot_1700 +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Including observations posterior to 1700**#

summary(c3_excl1700 <- ivreg(lead.polity ~ city_pop_share +
                               ln_oworld_1750 + ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + 
                               ln_elevation_1750 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + 
                               ln_tropical_1750 + ln_tropical_1800 + 
                               ln_tropical_1850 + ln_tropical_1900 + 
                               ln_rugged_1750 + ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share +
                               ln_wpot_1750 + ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)


#**Including observations posterior to 1800**#

summary(c3_excl1750 <- ivreg(lead.polity ~ city_pop_share + 
                               ln_oworld_1800 + 
                               ln_oworld_1850 + ln_oworld_1900 + 
                               ln_elevation_1800 + ln_elevation_1850 + ln_elevation_1900 + ln_tropical_1800 + 
                               ln_tropical_1850 + ln_tropical_1900 +                                ln_rugged_1800 + ln_rugged_1850 + 
                               ln_rugged_1900
                             + as.factor(country) + as.factor(year) - 1
                             | . - city_pop_share +
                               ln_wpot_1800 + ln_wpot_1850 + ln_wpot_1900, data=Potato), diagnostics=TRUE)

#**Final results with corrected SE**#

c3_noexcl <- cluster.robust.se(c3_noexcl, c3_noexcl$model$`as.factor(country)`)
c3_excl1100 <- cluster.robust.se(c3_excl1100, c3_excl1100$model$`as.factor(country)`)
c3_excl1200 <- cluster.robust.se(c3_excl1200, c3_excl1200$model$`as.factor(country)`)
c3_excl1300 <- cluster.robust.se(c3_excl1300, c3_excl1300$model$`as.factor(country)`)
c3_excl1400 <- cluster.robust.se(c3_excl1400, c3_excl1400$model$`as.factor(country)`)
c3_excl1500 <- cluster.robust.se(c3_excl1500, c3_excl1500$model$`as.factor(country)`)
c3_excl1600 <- cluster.robust.se(c3_excl1600, c3_excl1600$model$`as.factor(country)`)
c3_excl1700 <- cluster.robust.se(c3_excl1700, c3_excl1700$model$`as.factor(country)`)
c3_excl1750 <- cluster.robust.se(c3_excl1750, c3_excl1750$model$`as.factor(country)`)

#**Graphical results of the robustness check**#

model_Frame <- data.frame(Variable = c(rep('Model 1', 9), rep('Model 2', 9), rep('Model 3', 9)
),
Coefficient = c(c1_noexcl[1,1], c1_excl1100[1,1], c1_excl1200[1,1], c1_excl1300[1,1], c1_excl1400[1,1], c1_excl1500[1,1], c1_excl1600[1,1], c1_excl1700[1,1], c1_excl1750[1,1], c2_noexcl[1,1], c2_excl1100[1,1], c2_excl1200[1,1], c2_excl1300[1,1], c2_excl1400[1,1], c2_excl1500[1,1], c2_excl1600[1,1], c2_excl1700[1,1], c2_excl1750[1,1], c3_noexcl[1,1], c3_excl1100[1,1], c3_excl1200[1,1], c3_excl1300[1,1], c3_excl1400[1,1], c3_excl1500[1,1], c3_excl1600[1,1], c3_excl1700[1,1], c3_excl1750[1,1]
),
SE = c(c1_noexcl[1,2], c1_excl1100[1,2], c1_excl1200[1,2], c1_excl1300[1,2], c1_excl1400[1,2], c1_excl1500[1,2], c1_excl1600[1,2], c1_excl1700[1,2], c1_excl1750[1,2], c2_noexcl[1,2], c2_excl1100[1,2], c2_excl1200[1,2], c2_excl1300[1,2], c2_excl1400[1,2], c2_excl1500[1,2], c2_excl1600[1,2], c2_excl1700[1,2], c2_excl1750[1,2], c3_noexcl[1,2], c3_excl1100[1,2], c3_excl1200[1,2], c3_excl1300[1,2], c3_excl1400[1,2], c3_excl1500[1,2], c3_excl1600[1,2], c3_excl1700[1,2], c3_excl1750[1,2]
),
Period = c(' All', '>1100', '>1200', '>1300', '>1400', '>1500', '>1600', '>1700', '>1750', ' All', '>1100', '>1200', '>1300', '>1400', '>1500', '>1600', '>1700', '>1750', ' All', '>1100', '>1200', '>1300', '>1400', '>1500', '>1600', '>1700', '>1750'
),
show0 = c(rep(TRUE, 27)))

allModelFrame <- data.frame(rbind(model_Frame)) 

interval1 <- -qnorm((1-0.9))  # # Width of CI: 90% multiplier for one-tailed

allModelFrame$Variable <- factor(allModelFrame$Variable)

fig <- ggplot(allModelFrame, aes(group=Period, colour=Period, alpha = show0))
fig <- fig + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + labs(x="",y="Estimated Effects of Urban Population Share on \nDemocracy Score")
fig <- fig + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval1, ymax = Coefficient + SE*interval1),                              lwd = 0.75, position = position_dodge(width = 1/2),                              shape = 20)
fig <- fig + theme_bw() + theme(axis.text=element_text(size=14),
                                axis.title=element_text(size=18)) 
fig_appD <- fig + ggtitle("") + scale_alpha_identity()
print(fig_appD)

pdf("Figure.pdf", width = 10, height = 6)
fig_appD
dev.off()

#******************************#
#         Appendix E:          #
#      Sensitivity Analysis    #
#******************************#

#***Move to STATA .ado file***#




